'''
Created on Sep 17, 2010
@author: oabalbin

Cut the reads in pieces of 24 nucleotides
'''
from optparse import OptionParser
from collections import deque, defaultdict
 
class myInterval:
    def __init__(self,chr, start, end, name, score, strand, nexons=0, exonsizes=[], exonstart=[]):
        self.chr, self.start, self.end, self.name, self.score, self.strand =  \
        chr, int(start), int(end), name, score, strand

        if nexons:
            self.nexons, self.exonssize, self.exonstart = nexons, exonsizes, exonstart


def create_list_of_intervals(intervals_bed_file):
    
    intervals_bed_file = open(intervals_bed_file)
    interval_list=defaultdict()
    for line in intervals_bed_file:
        fields = line.strip('\n').split('\t')
        nfields = len(fields)
        
        if nfields < 3:
            print "The bed file has not enough information"
        elif nfields == 6:
            chr_tmp=fields[0].replace('chr','')
            name, score, strand=fields[3],fields[4],default_strand(fields[5])
            interval_list[name]=((myInterval(chr_tmp, int(fields[1]), int(fields[2]), name, score, strand)))
        elif nfields < 6:
            chr_tmp=fields[0].replace('chr','')
            name, score, strand="noname","noscore","+"
            interval_list[name]=((myInterval(chr_tmp, int(fields[1]), int(fields[2]), name, score, strand)))
        elif nfields == 12:
            chr_tmp=fields[0].replace('chr','')
            name, score, strand=fields[3],fields[4],default_strand(fields[5])
            nexons, exonsizes, exonstart = fields[9], fields[10].split(','),fields[11].split(',')
            interval_list[name]=((myInterval(chr_tmp, int(fields[1]), int(fields[2]), name, score, strand,\
                                             nexons, exonsizes, exonstart)))

    intervals_bed_file.close()
    return interval_list

def default_strand(x):
    if not x:
        return '+'
    elif x=='.':
        return '+'
    else:
        return x


def intersect_intervals_tree(list_of_intervals1, list_of_intervals2, outfile_name):
    
    outfile = open(outfile_name,'w')
     
    intersected_intervals=[]
    score=999
    
    for name in list_of_intervals1:        
        try:
            thisintv = list_of_intervals2[name]
            for i in range(int(thisintv.nexons)):
                estart = int(thisintv.start)+int(thisintv.exonstart[i])
                eend = estart + int(thisintv.exonssize[i])
                line = ["chr"+thisintv.chr,estart, eend, thisintv.name, score,thisintv.strand]   
                outfile.write(",".join(map(str,line)).replace(',','\t')+'\n')
        except KeyError:
            continue
    
    return intersected_intervals

if __name__ == '__main__':

    optionparser = OptionParser("usage: %prog [options] ")
    optionparser.add_option("-t", "--transcripts", dest="transcripts",
                            help="Bed file with the location of transcripts of interest")
    optionparser.add_option("-b", "--transexons", dest="transexons",
                            help="bed file with transcripts-exons of interest")
    
    optionparser.add_option("-o", "--exons", dest="exons",
                            help="output bed file with the exons of interest")

    (options, args) = optionparser.parse_args()

    
    interval_file1 = options.transcripts
    interval_file2 = options.transexons
    outputfile_name = options.exons
    interval_list1 = create_list_of_intervals(interval_file1)
    interval_list2 = create_list_of_intervals(interval_file2)
    intersect_intervals_tree(interval_list1, interval_list2, outputfile_name)
    
